# Load necessary libraries
library(forecast)
library(tseries)
library(TSA)
library(tidyverse)
library(rugarch)
# Load the data
aapl_monthly_data <- read.csv("~/Documents/GitHub/MA-641-Course-Project/AAPL_Monthly2016.csv")
# Convert the date column to Date type
aapl_monthly_data$Date <- as.Date(aapl_monthly_data$Date, format="%Y-%m-%d")
# Inspect the data
head(aapl_monthly_data)
summary(aapl_monthly_data)
Date Open High Low Close
Min. :2016-01-01 Min. : 23.49 Min. : 24.72 Min. : 22.37 Min. : 23.43
1st Qu.:2018-02-01 1st Qu.: 41.74 1st Qu.: 44.30 1st Qu.: 40.16 1st Qu.: 41.95
Median :2020-03-01 Median : 71.56 Median : 81.06 Median : 64.09 Median : 73.45
Mean :2020-03-01 Mean : 94.47 Mean :101.04 Mean : 88.97 Mean : 95.93
3rd Qu.:2022-04-01 3rd Qu.:148.99 3rd Qu.:157.50 3rd Qu.:138.27 3rd Qu.:149.80
Max. :2024-05-01 Max. :196.24 Max. :199.62 Max. :187.45 Max. :196.45
Adj.Close Volume
Min. : 21.39 Min. :9.697e+08
1st Qu.: 39.72 1st Qu.:1.676e+09
Median : 71.54 Median :2.240e+09
Mean : 93.98 Mean :2.320e+09
3rd Qu.:147.50 3rd Qu.:2.801e+09
Max. :195.41 Max. :6.280e+09
About the Data:
# Create a time series object
aaplmonthly_ts <- ts(aapl_monthly_data$Close, start=c(2016, 01), end = c(2024, 05), frequency=12)
# Descriptive Analysis
plot(aaplmonthly_ts, main="Monthly Apple Stock Prices", ylab="Close Price", xlab="Time")
summary(aaplmonthly_ts)
Min. 1st Qu. Median Mean 3rd Qu. Max.
23.43 41.95 73.45 95.93 149.80 196.45
boxplot(aaplmonthly_ts ~ cycle(aaplmonthly_ts), main="Seasonal Boxplot of Monthly Apple Stock Prices", ylab="Close Price")
Time Series Plot:
Summary:
Seasonal Boxplot:
# ACF and PACF Plots
par(mar=c(5, 5, 4, 2) + 0.1)
acf(aaplmonthly_ts, main="ACF of Monthly Apple Stock Prices", lag.max = 72)
pacf(aaplmonthly_ts, main="PACF of Monthly Apple Stock Prices", lag.max = 72)
eacf(aaplmonthly_ts)
AR/MA
0 1 2 3 4 5 6 7 8 9 10 11 12 13
0 x x x x x x x x x x x x x x
1 o x o o o o o o o o o o o o
2 x x o o o o o o x o o o o o
3 o o o o o o o o o o o o o o
4 x o o o o o o o o o o o o o
5 o o o o o o o o o o o o o o
6 o o o o o o o o o o o o o o
7 o o o x o o o o o o o o o o
ACF Plot:
PACF Plot:
# Augmented Dickey-Fuller Test
adf_test <- adf.test(aaplmonthly_ts, alternative="stationary")
print(adf_test)
Augmented Dickey-Fuller Test
data: aaplmonthly_ts
Dickey-Fuller = -2.3411, Lag order = 4, p-value = 0.4353
alternative hypothesis: stationary
# Differencing the series if it is not stationary
if (adf_test$p.value > 0.05) {
ts_data_diff <- diff(aaplmonthly_ts, differences=1)
adf_test_diff <- adf.test(ts_data_diff, alternative="stationary")
print(adf_test_diff)
# Update the time series data to the differenced series
aaplmonthly_ts <- ts_data_diff
}
Warning: p-value smaller than printed p-value
Augmented Dickey-Fuller Test
data: ts_data_diff
Dickey-Fuller = -5.0114, Lag order = 4, p-value = 0.01
alternative hypothesis: stationary
# Time Series Plot after Differencing
plot(aaplmonthly_ts, main="Monthly Apple Stock Prices", ylab="Close Price", xlab="Time")
# Fit AR model
ar_model <- Arima(aaplmonthly_ts, order=c(2,0,0))
summary(ar_model)
Series: aaplmonthly_ts
ARIMA(2,0,0) with non-zero mean
Coefficients:
ar1 ar2 mean
0.0409 -0.2708 1.6476
s.e. 0.0988 0.0980 0.6883
sigma^2 = 73.24: log likelihood = -355.13
AIC=718.27 AICc=718.69 BIC=728.69
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set -0.000308276 8.428628 6.23919 507.0341 563.7521 0.7036177 -0.001614169
# Perform diagnostics for AR(2)
tsdiag(ar_model, gof.lag = 10, main = "Diagnostics for AR(2)")
# Q-Q plot for AR(2)
residuals_ar2 <- residuals(ar_model)
qqnorm(residuals_ar2, main = "Q-Q Plot of Residuals for AR(2)")
qqline(residuals_ar2, col = "red")
# Fit MA model
ma_model <- Arima(aaplmonthly_ts, order=c(0,0,2))
summary(ma_model)
Series: aaplmonthly_ts
ARIMA(0,0,2) with non-zero mean
Coefficients:
ma1 ma2 mean
0.0231 -0.2775 1.6525
s.e. 0.0979 0.0971 0.6327
sigma^2 = 73.16: log likelihood = -355.09
AIC=718.17 AICc=718.59 BIC=728.59
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set -0.006781185 8.424279 6.24754 520.9368 585.7699 0.7045594 0.0114845
# Perform diagnostics for MA(2)
tsdiag(ma_model, gof.lag = 10, main = "Diagnostics for MA(2)")
# Q-Q plot for MA(2)
residuals_ma2 <- residuals(ma_model)
qqnorm(residuals_ma2, main = "Q-Q Plot of Residuals for MA(2)")
qqline(residuals_ma2, col = "red")
# Fit ARMA(1,1,1) model
arma_model1 <- Arima(aaplmonthly_ts, order=c(1,1,1))
summary(arma_model1)
Series: aaplmonthly_ts
ARIMA(1,1,1)
Coefficients:
ar1 ma1
0.0413 -1.0000
s.e. 0.1034 0.0308
sigma^2 = 78.94: log likelihood = -357.98
AIC=721.96 AICc=722.21 BIC=729.74
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.5179315 8.750795 6.505048 433.929 441.6441 0.7335995 -0.0008596899
# Perform diagnostics for ARIMA(1,1,1)
tsdiag(arma_model1, gof.lag = 10, main = "Diagnostics for ARIMA(1,1,1)")
# Q-Q plot for ARIMA(1,1,1)
residuals_arma111 <- residuals(arma_model1)
qqnorm(residuals_arma111, main = "Q-Q Plot of Residuals for ARIMA(1,1,1)")
qqline(residuals_arma111, col = "red")
# ARMA(2,1,1) Model
arma_model2 <- Arima(aaplmonthly_ts, order=c(2,1,1))
summary(arma_model2)
Series: aaplmonthly_ts
ARIMA(2,1,1)
Coefficients:
ar1 ar2 ma1
0.0486 -0.2632 -1.0000
s.e. 0.0996 0.0988 0.0377
sigma^2 = 74.02: log likelihood = -354.58
AIC=717.16 AICc=717.59 BIC=727.54
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.6253754 8.429532 6.225648 78.44247 127.6236 0.7020905 -0.01202426
# Perform diagnostics for ARIMA(2,1,1)
tsdiag(arma_model2, gof.lag = 10, main = "Diagnostics for ARIMA(2,1,1)")
# Q-Q plot for ARIMA(2,1,1)
residuals_arma211 <- residuals(arma_model2)
qqnorm(residuals_arma211, main = "Q-Q Plot of Residuals for ARIMA(2,1,1)")
qqline(residuals_arma211, col = "red")
# ARMA(2,1,2) Model
arma_model3 <- Arima(aaplmonthly_ts, order=c(2,1,2))
summary(arma_model3)
Series: aaplmonthly_ts
ARIMA(2,1,2)
Coefficients:
ar1 ar2 ma1 ma2
0.0375 -0.2628 -0.9880 -0.0120
s.e. 0.4098 0.1004 0.4284 0.4267
sigma^2 = 74.8: log likelihood = -354.58
AIC=719.16 AICc=719.8 BIC=732.13
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.6240689 8.429788 6.223491 79.15774 126.6565 0.7018472 -0.01277316
# Perform diagnostics for ARIMA(2,1,2)
tsdiag(arma_model3, gof.lag = 10, main = "Diagnostics for ARIMA(2,1,2)")
# Q-Q plot for ARIMA(2,1,2)
residuals_arma212 <- residuals(arma_model3)
qqnorm(residuals_arma212, main = "Q-Q Plot of Residuals for ARIMA(2,1,2)")
qqline(residuals_arma212, col = "red")
# Create a time series object
aaplmonthly_ts2 <- ts(aapl_monthly_data$Close, start=c(2016, 01), end = c(2024, 05), frequency=12)
# Calculate returns for modeling
returns <- diff(log(aaplmonthly_ts2))
returns <- returns[!is.na(returns)]
plot(returns, main="Monthly Apple Stock Prices", ylab="Close Price", xlab="Time", type="l")
# Display the fit summary
fit_garch
*---------------------------------*
* GARCH Model Fit *
*---------------------------------*
Conditional Variance Dynamics
-----------------------------------
GARCH Model : sGARCH(1,1)
Mean Model : ARFIMA(1,0,0)
Distribution : norm
Optimal Parameters
------------------------------------
Estimate Std. Error t value Pr(>|t|)
mu 0.020696 0.008198 2.52440 0.01159
ar1 0.009283 0.100840 0.09206 0.92665
omega 0.000212 0.000204 1.03801 0.29927
alpha1 0.000000 0.014482 0.00000 1.00000
beta1 0.969597 0.051349 18.88262 0.00000
Robust Standard Errors:
Estimate Std. Error t value Pr(>|t|)
mu 0.020696 0.006649 3.112732 0.001854
ar1 0.009283 0.103105 0.090038 0.928257
omega 0.000212 0.000379 0.558168 0.576730
alpha1 0.000000 0.005313 0.000000 1.000000
beta1 0.969597 0.059723 16.234816 0.000000
LogLikelihood : 108.877
Information Criteria
------------------------------------
Akaike -2.0775
Bayes -1.9473
Shibata -2.0822
Hannan-Quinn -2.0248
Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
statistic p-value
Lag[1] 0.0001857 0.9891
Lag[2*(p+q)+(p+q)-1][2] 1.2346095 0.5948
Lag[4*(p+q)+(p+q)-1][5] 2.4164966 0.5903
d.o.f=1
H0 : No serial correlation
Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
statistic p-value
Lag[1] 1.449 0.2287
Lag[2*(p+q)+(p+q)-1][5] 3.745 0.2876
Lag[4*(p+q)+(p+q)-1][9] 6.376 0.2572
d.o.f=2
Weighted ARCH LM Tests
------------------------------------
Statistic Shape Scale P-Value
ARCH Lag[3] 3.543 0.500 2.000 0.05978
ARCH Lag[5] 3.661 1.440 1.667 0.20717
ARCH Lag[7] 5.696 2.315 1.543 0.16300
Nyblom stability test
------------------------------------
Joint Statistic: 1.4525
Individual Statistics:
mu 0.07109
ar1 0.14752
omega 0.14307
alpha1 0.12991
beta1 0.14386
Asymptotic Critical Values (10% 5% 1%)
Joint Statistic: 1.28 1.47 1.88
Individual Statistic: 0.35 0.47 0.75
Sign Bias Test
------------------------------------
Adjusted Pearson Goodness-of-Fit Test:
------------------------------------
group statistic p-value(g-1)
1 20 20.4 0.3709
2 30 21.8 0.8284
3 40 34.4 0.6796
4 50 42.0 0.7504
Elapsed time : 0.05360389
# Specify GARCH(1,1) model
spec_garch <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1, 1)),
mean.model = list(armaOrder = c(2, 1)),
distribution.model = "norm")
fit_garch <- ugarchfit(spec = spec_garch, data = returns)
# Display the fit summary
fit_garch
*---------------------------------*
* GARCH Model Fit *
*---------------------------------*
Conditional Variance Dynamics
-----------------------------------
GARCH Model : sGARCH(1,1)
Mean Model : ARFIMA(2,0,1)
Distribution : norm
Optimal Parameters
------------------------------------
Estimate Std. Error t value Pr(>|t|)
mu 0.020739 0.005744 3.61058 0.000306
ar1 0.604229 0.341075 1.77154 0.076471
ar2 -0.151955 0.108878 -1.39564 0.162824
ma1 -0.610786 0.339011 -1.80167 0.071598
omega 0.000208 0.000238 0.87143 0.383518
alpha1 0.000000 0.013219 0.00000 1.000000
beta1 0.968842 0.047477 20.40660 0.000000
Robust Standard Errors:
Estimate Std. Error t value Pr(>|t|)
mu 0.020739 0.005612 3.69541 0.000220
ar1 0.604229 0.323349 1.86866 0.061671
ar2 -0.151955 0.095674 -1.58826 0.112228
ma1 -0.610786 0.293617 -2.08021 0.037506
omega 0.000208 0.000316 0.65791 0.510598
alpha1 0.000000 0.003947 0.00000 1.000000
beta1 0.968842 0.051773 18.71330 0.000000
LogLikelihood : 110.5264
Information Criteria
------------------------------------
Akaike -2.0705
Bayes -1.8882
Shibata -2.0795
Hannan-Quinn -1.9967
Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
statistic p-value
Lag[1] 0.0007236 0.9785
Lag[2*(p+q)+(p+q)-1][8] 0.7800872 1.0000
Lag[4*(p+q)+(p+q)-1][14] 1.9624525 0.9999
d.o.f=3
H0 : No serial correlation
Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
statistic p-value
Lag[1] 2.977 0.08444
Lag[2*(p+q)+(p+q)-1][5] 4.546 0.19337
Lag[4*(p+q)+(p+q)-1][9] 6.847 0.21233
d.o.f=2
Weighted ARCH LM Tests
------------------------------------
Statistic Shape Scale P-Value
ARCH Lag[3] 2.326 0.500 2.000 0.1273
ARCH Lag[5] 2.329 1.440 1.667 0.4030
ARCH Lag[7] 4.181 2.315 1.543 0.3212
Nyblom stability test
------------------------------------
Joint Statistic: 2.2869
Individual Statistics:
mu 0.14030
ar1 0.06721
ar2 0.04291
ma1 0.07467
omega 0.13026
alpha1 0.10397
beta1 0.13095
Asymptotic Critical Values (10% 5% 1%)
Joint Statistic: 1.69 1.9 2.35
Individual Statistic: 0.35 0.47 0.75
Sign Bias Test
------------------------------------
Adjusted Pearson Goodness-of-Fit Test:
------------------------------------
group statistic p-value(g-1)
1 20 27.6 0.091435
2 30 40.4 0.077606
3 40 47.2 0.172404
4 50 77.0 0.006497
Elapsed time : 0.0871799
# Plot diagnostics
plot(fit_garch)
Make a plot selection (or 0 to exit):
1: Series with 2 Conditional SD Superimposed 2: Series with 1% VaR Limits
3: Conditional SD (vs |returns|) 4: ACF of Observations
5: ACF of Squared Observations 6: ACF of Absolute Observations
7: Cross Correlation 8: Empirical Density of Standardized Residuals
9: QQ-Plot of Standardized Residuals 10: ACF of Standardized Residuals
11: ACF of Squared Standardized Residuals 12: News-Impact Curve
1
Make a plot selection (or 0 to exit):
1: Series with 2 Conditional SD Superimposed 2: Series with 1% VaR Limits
3: Conditional SD (vs |returns|) 4: ACF of Observations
5: ACF of Squared Observations 6: ACF of Absolute Observations
7: Cross Correlation 8: Empirical Density of Standardized Residuals
9: QQ-Plot of Standardized Residuals 10: ACF of Standardized Residuals
11: ACF of Squared Standardized Residuals 12: News-Impact Curve
2
please wait...calculating quantiles...
Make a plot selection (or 0 to exit):
1: Series with 2 Conditional SD Superimposed 2: Series with 1% VaR Limits
3: Conditional SD (vs |returns|) 4: ACF of Observations
5: ACF of Squared Observations 6: ACF of Absolute Observations
7: Cross Correlation 8: Empirical Density of Standardized Residuals
9: QQ-Plot of Standardized Residuals 10: ACF of Standardized Residuals
11: ACF of Squared Standardized Residuals 12: News-Impact Curve
3
Make a plot selection (or 0 to exit):
1: Series with 2 Conditional SD Superimposed 2: Series with 1% VaR Limits
3: Conditional SD (vs |returns|) 4: ACF of Observations
5: ACF of Squared Observations 6: ACF of Absolute Observations
7: Cross Correlation 8: Empirical Density of Standardized Residuals
9: QQ-Plot of Standardized Residuals 10: ACF of Standardized Residuals
11: ACF of Squared Standardized Residuals 12: News-Impact Curve
4
Make a plot selection (or 0 to exit):
1: Series with 2 Conditional SD Superimposed 2: Series with 1% VaR Limits
3: Conditional SD (vs |returns|) 4: ACF of Observations
5: ACF of Squared Observations 6: ACF of Absolute Observations
7: Cross Correlation 8: Empirical Density of Standardized Residuals
9: QQ-Plot of Standardized Residuals 10: ACF of Standardized Residuals
11: ACF of Squared Standardized Residuals 12: News-Impact Curve
5
Make a plot selection (or 0 to exit):
1: Series with 2 Conditional SD Superimposed 2: Series with 1% VaR Limits
3: Conditional SD (vs |returns|) 4: ACF of Observations
5: ACF of Squared Observations 6: ACF of Absolute Observations
7: Cross Correlation 8: Empirical Density of Standardized Residuals
9: QQ-Plot of Standardized Residuals 10: ACF of Standardized Residuals
11: ACF of Squared Standardized Residuals 12: News-Impact Curve
6
Make a plot selection (or 0 to exit):
1: Series with 2 Conditional SD Superimposed 2: Series with 1% VaR Limits
3: Conditional SD (vs |returns|) 4: ACF of Observations
5: ACF of Squared Observations 6: ACF of Absolute Observations
7: Cross Correlation 8: Empirical Density of Standardized Residuals
9: QQ-Plot of Standardized Residuals 10: ACF of Standardized Residuals
11: ACF of Squared Standardized Residuals 12: News-Impact Curve
7
Make a plot selection (or 0 to exit):
1: Series with 2 Conditional SD Superimposed 2: Series with 1% VaR Limits
3: Conditional SD (vs |returns|) 4: ACF of Observations
5: ACF of Squared Observations 6: ACF of Absolute Observations
7: Cross Correlation 8: Empirical Density of Standardized Residuals
9: QQ-Plot of Standardized Residuals 10: ACF of Standardized Residuals
11: ACF of Squared Standardized Residuals 12: News-Impact Curve
8
Make a plot selection (or 0 to exit):
1: Series with 2 Conditional SD Superimposed 2: Series with 1% VaR Limits
3: Conditional SD (vs |returns|) 4: ACF of Observations
5: ACF of Squared Observations 6: ACF of Absolute Observations
7: Cross Correlation 8: Empirical Density of Standardized Residuals
9: QQ-Plot of Standardized Residuals 10: ACF of Standardized Residuals
11: ACF of Squared Standardized Residuals 12: News-Impact Curve
9
Make a plot selection (or 0 to exit):
1: Series with 2 Conditional SD Superimposed 2: Series with 1% VaR Limits
3: Conditional SD (vs |returns|) 4: ACF of Observations
5: ACF of Squared Observations 6: ACF of Absolute Observations
7: Cross Correlation 8: Empirical Density of Standardized Residuals
9: QQ-Plot of Standardized Residuals 10: ACF of Standardized Residuals
11: ACF of Squared Standardized Residuals 12: News-Impact Curve
10
Make a plot selection (or 0 to exit):
1: Series with 2 Conditional SD Superimposed 2: Series with 1% VaR Limits
3: Conditional SD (vs |returns|) 4: ACF of Observations
5: ACF of Squared Observations 6: ACF of Absolute Observations
7: Cross Correlation 8: Empirical Density of Standardized Residuals
9: QQ-Plot of Standardized Residuals 10: ACF of Standardized Residuals
11: ACF of Squared Standardized Residuals 12: News-Impact Curve
11
Make a plot selection (or 0 to exit):
1: Series with 2 Conditional SD Superimposed 2: Series with 1% VaR Limits
3: Conditional SD (vs |returns|) 4: ACF of Observations
5: ACF of Squared Observations 6: ACF of Absolute Observations
7: Cross Correlation 8: Empirical Density of Standardized Residuals
9: QQ-Plot of Standardized Residuals 10: ACF of Standardized Residuals
11: ACF of Squared Standardized Residuals 12: News-Impact Curve
12
Make a plot selection (or 0 to exit):
1: Series with 2 Conditional SD Superimposed 2: Series with 1% VaR Limits
3: Conditional SD (vs |returns|) 4: ACF of Observations
5: ACF of Squared Observations 6: ACF of Absolute Observations
7: Cross Correlation 8: Empirical Density of Standardized Residuals
9: QQ-Plot of Standardized Residuals 10: ACF of Standardized Residuals
11: ACF of Squared Standardized Residuals 12: News-Impact Curve
0
# Forecasting with the TGARCH model
forecast_garch <- ugarchforecast(fit_garch, n.ahead=12)
plot(forecast_garch, which=1) # Forecast series
# Comparing Models using AIC and BIC
models <- list(ar_model, ma_model, arma_model1, arma_model2, arma_model3)
model_names <- c("AR", "MA", "ARIMA(1,1,1)", "ARIMA(2,1,1)", "ARIMA(2,1,2", "GARCH")
aic_values <- c(sapply(models, AIC), -2.0775)
bic_values <- c(sapply(models, BIC), -1.9473)
comparison <- data.frame(Model=model_names, AIC=aic_values, BIC=bic_values)
print(comparison)